/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*  Copyleft (!c) 2004 Div. Nucl. Med., Cyric, Tohoku Univ.         *
*  Copyleft (!c) 2004-2011 Elseif Laboratory.                      *
*  elseifkk@users.sf.net                                           *
*                                                                  *
*  All Rights Reversed.                                            *
*                                                                  *
*                                                                  *
*  LICENSE AGREEMENT:                                              *
*                                                                  *
*  You agree that we don't agree to anything.                      *
*                                                                  *
*                                                                  *
*  DISCLAIMER OF WARRANTIES:                                       *
*                                                                  *
*  This software is provided to you "AS IS," and the authors and   *
*  contributors disclaim any warranty or liability obligations to  *
*  you of any kind, implied or express, including but not limited  *
*  to the warranties of merchantability, fitness for a particular  *
*  purpose, and non-infringement.                                  *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <time.h>
#include "fbp_cti.h"
//---------------------------------------------
unsigned int  static _nX_ = 128; // image dimension in X
unsigned int  static _nY_ = 128; // image dimension in Y
unsigned int  static _jX0_;      // index of image origin in X
unsigned int  static _jY0_;      // index of image origin in Y
//-------------------------------------------------------------
double static y_i[_nt_][_nx_];
double static lmd_j[_nYmax_][_nXmax_];
short static ilmd[_nYmax_][_nXmax_];
double static fq[_nt_][_nx_*_BFS_];
double static fx[_nfx_];
char static headers[3][_BS_];
double static Delta; // sampling width
double static g_width;
double static win_slope;
//-------------------------------------------------------------
// header info and all
short static num_planes;
short static num_frames;
short static num_gates;
short static num_bed_position;
int static total_number_of_images;
int static i_bed_position;
int static i_planes;
int static i_frames;
int static i_phase;
int static i_record;
int static i_matrix;
int static last_record=0;
//
double static sample_distance;
double static pixel_size;
double static slice_width;
short static image_max,image_min;
float static quant_scale;
float static recon_scale;
float static isotope_halflife;
int static frame_start_time;
int static frame_duration;
float static loss_correction_factor;
short static dimension_1;
short static dimension_2;
short static filter_code;
float static decay_correction_factor;
//
int static number_of_matrix_created;
int static number_of_matrix_available;
int static forward_block;
int static backward_block;;
//
double static scale_factor;
char static original_file_name[32];
int static recon_duration;
time_t static t0;
//---------------------------------------------
void init_const()
{
  _jX0_      = _nX_/2;
  _jY0_      = _nY_/2;
  return;
}
//---------------------------------------------
void lower(long p, int l)
{
  int i;
  char *cp;
  for(i=0;i<l;i++){
    cp=(char *)(p+i);
    if(*cp==0) break;
    *cp=*cp|0x20;
  }
  return;
}
//---------------------------------------------
void f2cstr(long p, int l)
{
  int i;
  char *cp;
  for(i=0;i<l;i++){
    cp=(char *)(p+i);
    if(*cp==0x20){
      *cp=0;
    }
  }
  return;
}
//---------------------------------------------
int at_exp(long p, int l)
{
  int i;
  char *cp;
  for(i=l-1;0<=i;i--){
    cp=(char *)(p+i);
    if(*cp==0x2E) return i;
  }
  return 0;
}
//---------------------------------------------
void clelmd()
{
  unsigned int jX,jY;
  for(jY=0;jY<_nY_;jY++){
    for(jX=0;jX<_nX_;jX++){
      lmd_j[jY][jX]=0;
    }
  }
  return;
}
//---------------------------------------------
double gauss(double x)
{
  return exp(-pow(x/g_width,2.0)/2.0)/(g_width*sqrt(PI2));
}
//---------------------------------------------
void smtyi()
{
  int const N=4;
  double h[N+2];
  double s;
  int j;
  int ix,it;
  double x;
  double norm;
  double y1,y2;
  double yy[_nx_];
  norm=0.0;
  for(j=0;j<=N;j++){
    x=(double)(j)*Delta;
    h[j]=gauss(x);
    norm=norm+h[j];
  }
  norm=1.0/(2.0*norm-h[0]);
  for(it=0;it<_nt_;it++){
    for(ix=0;ix<_nx_;ix++){
      s=0.0;
      for(j=N;1<=j;j--){
	if(ix-j<0){
	  y1=0.0;
	}else{
	  y1=y_i[it][ix-j];
	}
	if(_nx_<ix+j){
	  y2=0.0;
	}else{
	  y2=y_i[it][ix+j];
	}
	s=s+h[j]*(y1+y2);
      }
      yy[ix]=(s+h[0]*y_i[it][ix])*norm;
    }
    for(ix=0;ix<_nx_;ix++){
      y_i[it][ix]=yy[ix];
    }
  }
  return;
}
//---------------------------------------------
int scn2mtv(FILE *f, unsigned int nx, unsigned int nt, int off1, int off2, int n)
{
  unsigned int offset;
  short w;
  unsigned int ix,it;
  offset = off1+(n-1)*((nx*nt*2)+off2)+(i_record-1)*_BS_;
  if(fseek(f,offset-_BS_,SEEK_SET)!=0) return 1;
  if(fread(&headers[HID_SUB][0],_BS_,1,f)!=1) return 2;
  if(fseek(f,offset,SEEK_SET)!=0) return 3;
  for(it=0;it<nt;it++){
    for(ix=0;ix<nx;ix++){
      if(fread(&w,2,1,f)!=1) return 3;
      if(w<0) w=0;
      y_i[it][ix]=(double)w;
    }
  }
  return 0;
}
//---------------------------------------------
double fold_Ram_Lak(int m)
{
  if(m==0) return PI/(2.0*pow(Delta,2.0));
  if(m & 1) return -2.0/(PI*pow((double)m*Delta,2.0));
  return 0.0;
}
//---------------------------------------------
double fold_Ram_Lak_Win(int m)
{
  double a;
  if(m==-1||m==1){
    a=pow(PIH/Delta,2.0);
  }else if(m & 1){ // odd
    a=0.0;
  }else{
    a=-2.0*(1.0+pow((double)m,2.0))/pow(Delta*(1.0-pow((double)m,2.0)),2.0);
  }
  return fold_Ram_Lak(m)*win_slope+(1.0-win_slope)*a/(PI);
}
//---------------------------------------------
double fold_Shepp_Logan(int m)
{
  return 4.0/(PI*pow(Delta,2.0)*(1.0-4.0*pow((double)m,2.0)));
}
//---------------------------------------------
double fold_Shepp_Logan_Win(int m)
{
  double a;
  a=-8.0*pow(1.0/Delta,2.0)*(1.0/(1.0-pow((double)(2*m),2.0))-3.0/(9.0-pow((double)(2*m),2.0)));
  return fold_Shepp_Logan(m)*win_slope+(1.0-win_slope)*a/(PI2);
}
//---------------------------------------------
double fold_None(int m)
{
  if(m==0) return 1.0/Delta;
  return 0.0;
}
//-----------------------------------------------------
void setfx(int idf)
{
  int m;
  switch(idf){
  case IDF_SHEPP_LOGAN:
    for(m=0;m<_nfx_;m++){
      fx[m]=fold_Shepp_Logan(m);
    }
    break;
  case IDF_RAM_LAK:
    for(m=0;m<_nfx_;m++){
      fx[m]=fold_Ram_Lak(m);
    }
    break;
  case IDF_NONE:
    for(m=0;m<_nfx_;m++){
      fx[m]=fold_None(m);
    }
    break;
  default:
    switch(idf){
    case IDF_SHEPP_LOGAN_HANNING:
    case IDF_RAM_LAK_HANNING:
      win_slope=0.54;
      break;
    case IDF_RAM_LAK_HAMMING:
    case IDF_SHEPP_LOGAN_HAMMING:
      win_slope=0.50;
      break;
    }
    if(idf==IDF_RAM_LAK_HAMMING||idf==IDF_RAM_LAK_HANNING){
      for(m=0;m<_nfx_;m++){
        fx[m]=fold_Ram_Lak_Win(m);
      }
    }else{
      for(m=0;m<_nfx_;m++){
        fx[m]=fold_Shepp_Logan_Win(m);
      }
    }
    switch(idf){
    case IDF_RAM_LAK_HANNING:
      win_slope=0.54;
      break;
    case IDF_RAM_LAK_HAMMING:
      win_slope=0.50;
      break;
    }
    for(m=0;m<_nfx_;m++){
      fx[m]=fold_Ram_Lak_Win(m);
    }
    break;
    switch(idf){
    case IDF_SHEPP_LOGAN_HANNING:
    case IDF_RAM_LAK_HANNING:
      win_slope=0.54;
      break;
    case IDF_RAM_LAK_HAMMING:
    case IDF_SHEPP_LOGAN_HAMMING:
      win_slope=0.50;
      break;
    }
    if(idf==IDF_RAM_LAK_HAMMING||idf==IDF_RAM_LAK_HANNING){
      for(m=0;m<_nfx_;m++){
        fx[m]=fold_Ram_Lak_Win(m);
      }
    }else{
      for(m=0;m<_nfx_;m++){
        fx[m]=fold_Shepp_Logan_Win(m);
      }
    }
    break;
  }  
  return;
}
//-----------------------------------------------------
void setfq()
{
  unsigned int ix,it;
  int m,mm;
  double s;
  for(ix=0;ix<_nx_*_BFS_;ix++){
    for(it=0;it<_nt_;it++){
      s=0.0;
      mm=(int)(ix-_ix0_*_BFS_);
      for(m=-_ix0_;m<=_ix0_;m++){
	s=s+fx[abs(mm-m)]*y_i[it][m+_ix0_];
      }
      fq[it][ix]=s;
    }
  }
  return;
}
//---------------------------------------------
void read_header(int hid, FILE *f)
{
  int off;
  switch(hid){
  case HID_MAIN:
    off=0;
    break;
  case HID_MATREC:
    off=_BS_;
    break;
  case HID_SUB:
    off=_BS_*2;
    break;
  default:
    return;
  }
  fseek(f,off,SEEK_SET);
  fread(&headers[hid][0],_BS_,1,f);
  return;
}
//---------------------------------------------
void write_header(int hid, FILE *f)
{
  int off;
  switch(hid){
  case HID_MAIN:
    off=0;
    break;
  case HID_MATREC:
    off=_BS_;
    break;
  case HID_SUB:
    off=_BS_*2;
    break;
  default:
    return;
  }
  fseek(f,off,SEEK_SET);
  fwrite(&headers[hid][0],_BS_,1,f);
  return;
}
//---------------------------------------------
void clear_header(int hid)
{
  int i;
  for(i=0;i<_BS_;i++){
    headers[hid][i]=0;
  }
  return;
}
//---------------------------------------------
double f2ilmd(double *min, double *max, short *wmin, short *wmax, int opt)
{
  unsigned int jX,jY;
  unsigned short const WMAX=0x7FFF; 
  double a;
  short w;
  *min=1.0e38;
  *max=-1.0e38;
  for(jY=0;jY<_nY_;jY++){
    for(jX=0;jX<_nX_;jX++){
      if(lmd_j[jY][jX]<*min) *min=lmd_j[jY][jX];
      if(*max<lmd_j[jY][jX]) *max=lmd_j[jY][jX];
    }
  }
  a=(double)WMAX/(*max);
  *wmin=WMAX;
  *wmax=-WMAX;
  for(jY=0;jY<_nY_;jY++){
    for(jX=0;jX<_nX_;jX++){
      w=(short)rint(lmd_j[jY][jX]*a);
      if(opt&&w<0) w=0;
      if(w<*wmin) *wmin=w;
      if(*wmax<w) *wmax=w;
      ilmd[jY][jX]=w;
    }
  }
  return 1.0/a;
}
//---------------------------------------------
void write_data_seq(FILE *f,int i)
{
  unsigned int jY;
  fwrite(&headers[HID_SUB][0],_BS_,1,f);
  for(jY=0;jY<_nY_;jY++){
    fwrite(&ilmd[jY][0],_DATABYTE_,_nX_,f);
  }
  return;
}
//---------------------------------------------
double v2f(short lo, short up)
{
  float *fp;
  int dw;
  fp=(float *)&dw;
  dw=((0x00FFFF&up)|(0xFFFF0000&(lo<<16)));
  return *fp/4.0;
}
//---------------------------------------------
void f2v(float f, short *up, short *lo)
{
  float *fp;
  int dw;
  fp=(float *)&dw;
  *fp=f*4.0;
  *lo=0x0000FFFF&dw;
  *up=(0xFFFF0000&dw)>>16;
  return;
}
//---------------------------------------------
void fill_header_param_white(int n, int i, int len)
{
  char *cp;
  int j;
  for(j=0;j<len;j++){
    cp=&headers[n][i+j];
    *cp=(char)(0x20);
  }
  return;
}
//---------------------------------------------
void set_header_param_s(int n, int i, short x)
{
  char *cp;
  cp=&headers[n][i];
  *cp=(char)(x&0x0FF);
  cp=&headers[n][i+1];
  *cp=(char)((x>>8)&0x00FF);
  return;
}
//---------------------------------------------
void set_header_param_i(int n, int i, int x)
{
  char *cp;
  cp=&headers[n][i];
  *cp=(char)(x&0x0FF);
  cp=&headers[n][i+1];
  *cp=(char)((x>>8)&0x00FF);
  cp=&headers[n][i+2];
  *cp=(char)((x>>16)&0x00FF);
  cp=&headers[n][i+3];
  *cp=(char)((x>>24)&0x00FF);
  return;
}
//---------------------------------------------
void set_header_param_f(int n, int i, float x)
{
  short lo,up;
  char *cp;
  f2v(x,&lo,&up);
  cp=&headers[n][i];
  *cp=(char)(lo&0x00FF);
  cp=&headers[n][i+1];
  *cp=(char)((lo>>8)&0x00FF);
  cp=&headers[n][i+2];
  *cp=(char)(up&0x00FF);
  cp=&headers[n][i+3];
  *cp=(char)((up>>8)&0x00FF);
  return;
}
//---------------------------------------------
float get_header_param_f(int n, int i)
{
  short lo,up;
  lo=(headers[n][i]&0x00FF)|((headers[n][i+1]<<8)&0xFF00);
  up=(headers[n][i+2]&0x00FF)|((headers[n][i+3]<<8)&0xFF00);
  return v2f(lo,up);
}
//---------------------------------------------
void get_header_param_i(int n, int i, int *x)
{
  *x=(headers[n][i]&0x00FF)
    |((headers[n][i+1]<<8)&0xFF00)
    |((headers[n][i+2]<<16)&0xFF0000)
    |((headers[n][i+3]<<24)&0xFF000000);
  return;
}
//---------------------------------------------
void get_header_param_s(int n, int i, short *x)
{
  *x=(headers[n][i]&0x00FF)|((headers[n][i+1]<<8)&0xFF00);
  return;
}
//---------------------------------------------
void set_img_header(int hid)
{
  int m;
  int l;
  int i;
  long p;
  int base;
  struct tm *tm0;
  static int fill[] = {
    0x00010100,0xffff0016,0x003e0019,0x006f004f,
    0x000000a3, 0,0,0,
    0,0,0,0,
    0xffff0025,0,0,0,
    0,0,0,0,
    0,0,0,0,
    0,0,0,0,
    0,0,0x00110000,0x0000ffff,
  };
  switch(hid){
  case HID_MAIN:
    fill_header_param_white(hid,28,20);
    memcpy(&headers[hid][28],&original_file_name[0],strlen(original_file_name));
    headers[hid][54]=2; // FType_image 
    break;
  case HID_MATREC:
    base=forward_block+1;
    backward_block=forward_block;
    forward_block=forward_block+(1+_nX_*_nY_*2/_BS_)*number_of_matrix_created+1;
    if(last_record){
      set_header_param_i(hid,4,2);
    }else{
      set_header_param_i(hid,4,forward_block);
    }
    if(i_record==1){
      set_header_param_i(hid,8,0);
    }else{
      set_header_param_i(hid,8,backward_block);
    }
    for(i=1;i<=number_of_matrix_created;i++){
      set_header_param_i(hid,i*16+4,base);
      set_header_param_i(hid,i*16+8,base+(_nX_*_nY_*2/_BS_));
      base=base+(_nX_*_nY_*2/_BS_)+1;
    }
    break;
  case HID_SUB:
    clear_header(HID_SUB);
    memcpy(&headers[hid][0],fill,0x7e);
    headers[hid][126]=2;
    headers[hid][128]=2;
    set_header_param_s(hid,132,_nX_);
    set_header_param_s(hid,134,_nY_);
    headers[hid][156]=0x20;
    headers[hid][158]=0xff;
    headers[hid][159]=0xff;
    set_header_param_f(hid,168,recon_scale);
    set_header_param_f(hid,172,quant_scale);
    set_header_param_i(hid,176,image_min);
    set_header_param_i(hid,178,image_max);
    set_header_param_f(hid,184,pixel_size);
    set_header_param_f(hid,188,slice_width);
    set_header_param_i(hid,192,frame_duration);
    set_header_param_i(hid,196,frame_start_time);
    set_header_param_s(hid,200,0xff98);
    tm0=localtime(&t0);
    set_header_param_s(hid,202,(short)(tm0->tm_hour));
    set_header_param_s(hid,204,(short)(tm0->tm_min));
    set_header_param_s(hid,206,(short)(tm0->tm_sec));
    set_header_param_i(hid,208,recon_duration*1000);
    headers[hid][220]=0x34; // ?
    headers[hid][222]=0xff; // ?
    headers[hid][223]=0xff; // ?
    set_header_param_s(hid,236,filter_code);
    m=0x01000000|(0x00FF0000&(i_planes<<16))|(0xFF00&(i_bed_position<<12))|1;
    set_header_param_i(hid,238,m);
    set_header_param_i(hid,242,m&0xFFFF00FF);
    set_header_param_i(hid,246,m);
    set_header_param_i(hid,300,scale_factor);
    set_header_param_i(hid,300,0xd5d04092);
    set_header_param_f(hid,304,decay_correction_factor);
    set_header_param_f(hid,308,loss_correction_factor);
    set_header_param_i(hid,324,0xffff0020);
    headers[hid][380]=2; // ?
    set_header_param_s(hid,382,(short)(tm0->tm_mday));
    set_header_param_s(hid,384,(short)(tm0->tm_mon));
    set_header_param_s(hid,386,(short)(tm0->tm_year+1900));
    set_header_param_f(hid,388,1.0);
    set_header_param_f(hid,392,1.0);
    headers[hid][397]=0x40;
    headers[hid][405]=0x40;
    fill_header_param_white(hid,420,40);
    p=(long)&original_file_name[0];
    l=at_exp(p,strlen(original_file_name));
    if(l==0) l=strlen(original_file_name);
    if(20<l) l=20;
    memcpy(&headers[hid][420],&original_file_name[0],l);
    memcpy(&headers[hid][420+l+1],"FR",2);
    headers[hid][420+l+1+3]=0x30|i_frames;
    memcpy(&headers[hid][420+l+6],"PL",2);
    headers[hid][420+l+1+8]=0x30|i_planes;
    memcpy(&headers[hid][420+l+11],"BD",2);
    headers[hid][420+l+1+13]=0x30|i_bed_position;
    break;
  }
  return;
}
//---------------------------------------------
void echo_header_params(int hid)
{
  switch(hid){
  case HID_MAIN:
    fprintf(stderr,"Number of Planes = %d\n",num_planes);
    fprintf(stderr,"Number of Frames = %d\n",num_frames);
    fprintf(stderr,"Number of Gates = %d\n",num_gates);
    fprintf(stderr,"Number of Bed Positions = %d\n",num_bed_position);
    break;
  case HID_MATREC: 
    fprintf(stderr,"Number of available entries in this record = %d\n",number_of_matrix_available);
    fprintf(stderr,"Number of created entries in this record = %d\n",number_of_matrix_created);
    break;
  case HID_SUB:
    fprintf(stderr,"Number of Total Views = %d\n",dimension_1);
    fprintf(stderr,"Number of Total Elements = %d\n",dimension_2);
    fprintf(stderr,"Sample distance = %g cm\n",sample_distance);
    fprintf(stderr,"Isotope halflife = %g sec\n",isotope_halflife);
    fprintf(stderr,"Frame start time = %d msec\n",frame_start_time);
    fprintf(stderr,"Frame duration = %d msec\n",frame_duration);
    fprintf(stderr,"Scale factor = %g\n",scale_factor);
    fprintf(stderr,"Loss correction factor = %g\n",loss_correction_factor);
    break;
  }
  return;
}
//---------------------------------------------
void read_someheader_params(int ient)
{
  int base;
  base=ient*16;
  get_header_param_i(HID_MATREC,base,&i_frames);
  i_frames=i_frames&0xFFF;
  i_bed_position=headers[HID_MATREC][base+1]>>4;
  i_planes=headers[HID_MATREC][base+2];
  i_phase=headers[HID_MATREC][base+3]&0x3F;
  return;
}
//---------------------------------------------
void read_header_params(int hid)
{
  switch(hid){
  case HID_MAIN:
    memcpy(&original_file_name[0],&headers[hid][28],20);
    get_header_param_s(hid,376,&num_planes);
    get_header_param_s(hid,378,&num_frames);
    get_header_param_s(hid,380,&num_gates);
    get_header_param_s(hid,382,&num_bed_position);
    break;
  case HID_MATREC:
    get_header_param_i(hid,0,&number_of_matrix_available);
    get_header_param_i(hid,4,&last_record);
    if(last_record!=2) last_record=0;
    get_header_param_i(hid,12,&number_of_matrix_created);
    break;
  case HID_SUB:
    get_header_param_s(hid,132,&dimension_1);
    get_header_param_s(hid,134,&dimension_2);
    sample_distance=get_header_param_f(hid,146);
    isotope_halflife=get_header_param_f(hid,166);
    scale_factor=get_header_param_f(hid,182);
    get_header_param_i(hid,456,&frame_start_time);
    get_header_param_i(hid,460,&frame_duration);
    loss_correction_factor=get_header_param_f(hid,464);
    break;
  }
  return;
}
//---------------------------------------------
float get_decay_correction_factor()
{
  float tau;
  tau=isotope_halflife/log(2.0)*1000.0;
  return (float)frame_duration/tau
         *exp((float)frame_start_time/tau)
         /(1.0-exp(-(float)frame_duration/tau));
}
//---------------------------------------------
int guess_filesize(fid)
{
  int n;
  int m;
  n=_BS_*2;
  if(31<total_number_of_images){
    m=total_number_of_images/31*_BS_;
  }else{
    m=0;
  }
  switch(fid){
  case FID_IMAGE:
    return n=n+(_nX_*_nY_*2+_BS_)*total_number_of_images+m;
  case FID_SCAN:
    return n=n+(_nx_*_nt_*2+_BS_)*total_number_of_images+m;
  }
  return 0;
}
//---------------------------------------------
int guess_number_of_images()
{
  int n;
  n=1;
  if(num_planes!=0) n=n*num_planes;
  if(num_frames!=0) n=n*num_frames;
  if(num_gates!=0) n=n*num_gates;
  if(num_bed_position!=0) n=n*num_bed_position;
  return n;
}
//---------------------------------------------
void print_info()
{
  int i=0;
  static char *info[] = {
  " ",
  PROGRAM,
  " ",
  "Description:",
  "Creates PET images from attenuation and normalization corrected PET", 
  "sinograms of 192x256 pixels using the folded back-projection method.",
  " ",
  "Usage:",
  "fbp_cti [Options] Sinogram Output_image",
  " ",
  "Options:",  
  " * Filter code (default Hann):",
  "        -F{None|Ramp|Butter|Hann|Hamm|Parzen|Shepp|SHann|SHamm}",
  " * Zoom factor (default 1.0):",
  "        -Z{zoom}",
  " * Image rotation (positive for clockwise, default 0):",
  "        -R{rotation in degrees}",
  " * Image shift in x and y direction (default 0):",
  "        -X{shift in cm}, -Y{shift in cm},",
  "        -x{shift in pixel}, -y{shift in pixel}",
  " * Image dimension in pixels (it must be even and <=512,",
  "   default 128):",
  "        -D{dim}",
  " * Gaussian smoothing width (no smoothing by default):",
  "        -S{width in cm}, or -s{width in pixel}",
  " * Extra normalization factor; the results will be multiplied by" ,
  "   10^norm:",
  "        -N{norm}",
  " * etc:",
  " \t-I\t Use 4-point mean value for pixel value.",
  " \t-a\t Disable Scale_factor.",  
  " \t-c\t Keep case of the file name (usually in capital) when set",
  " \t  \t by -e option; otherwise lower case letter is used.",
  " \t-e\t Set Output_image automatically to Original_file_name of",
  " \t  \t main header info.",
  " \t-f\t Force overwrite if file already exists.",
  " \t-k\t Use absolute unit of cps/cc, otherwise the value will be",
  " \t  \t normalized by zoom and dimension.",
  " \t-n\t Allow Image_min of subheader info to be negative,",
  " \t  \t otherwise 0.",
  " \t-p\t Replace negative value with 0.",
  " \t-q\t Quiet mode.",
  " \t-u\t Disable decay correction.",
  " \t-v\t Verbose mode.",
  " ",
  0};
  while(info[i]!=0) printf("%s\n", info[i++]);
  return;
}
//---------------------------------------------
int main(int argc, char *argv[])
{
  /* Ground control */
  unsigned int mana=0;
  int idf=IDF_RAM_LAK_HANNING;
  double zoom=1.0;
  double rot=0.0;
  double DX=0.0,DY=0.0;
  double iDX=0.0,iDY=0.0;
  char *filterName[]={
  "None","Ramp","Butter","Hann","Hamm","Parzen","Shepp",
  "SHann","SHamm", 0
  };
  double Gods_will=1.0/16384.0; // extra normalization factor
  FILE *fin,*fout;
  char fnin[128],fnout[128];
  /* las ninas */
  unsigned int emi   = 0;
  unsigned int tra   = 0;
  unsigned int tyomo = 0;
  /* USS longshot */
  double const hm = (1.0-1.0/sqrt(3.0))/2.0;
  double const hp = (1.0+1.0/sqrt(3.0))/2.0;
  int fsin;
  int fsout;
  double dx[4],dy[4];
  int i2;
  double gw=1;
  int first=1;
  int n_in,n_out;
  /* misc */
  double s;
  int n;
  unsigned int jX,jY;
  unsigned int it;
  double a;
  double d,x;
  double kx,ky;
  int m0;
  double f0;
  double f0sqr;
  double r;
  double norm;
  double min,max;
  int i;
  char *pc;
  int mm;
  long p;
  time_t t1;
  double x0,y0;
  char str_total[8];

  /* Ground control */
  mana = 0;
  if(argc<2){
    fprintf(stderr,"*** Error: too few arguments");
    print_info();
    exit(1);
  }
  fnin[0]=fnout[0]=(char)0;
  for(i=1;i<argc;i++){
    if(*argv[i]=='-' && strlen(argv[i])>1){
      pc=argv[i]+2;
      if(*(argv[i]+1)=='Z'){
	zoom=atof(pc);
	mana = mana | ZOOM_SET;
	continue;
      }else if(*(argv[i]+1)=='R'){
	rot=atof(pc);
	mana = mana | ROT_SET;
	continue;
      }else if(*(argv[i]+1)=='X'){
	DX=atof(pc);
	mana = mana | DX_SET;
	continue;
      }else if(*(argv[i]+1)=='Y'){
	DY=atof(pc);
	mana = mana | DY_SET;
	continue;
      }else if(*(argv[i]+1)=='x'){
	iDX=atof(pc);
	mana = mana | IDX_SET;
	continue;
      }else if(*(argv[i]+1)=='y'){
	iDY=atof(pc);
	mana = mana | IDY_SET;
	continue;
      }else if(*(argv[i]+1)=='s'){
	gw=atof(pc);
	mana = mana & (~GWcm_SET);
	mana = mana | GWpx_SET;;
	if(gw==0.0){
	  fprintf(stderr,"*** Error: Gaussian width must not be 0.\n");
	  exit(1);
	}
	continue;
      }else if(*(argv[i]+1)=='S'){
	gw=atof(pc);
	mana = mana & (~GWpx_SET);
	mana = mana | GWcm_SET;
	if(gw==0.0){
	  fprintf(stderr,"*** Error: Gaussian width must not be 0.\n");
	  exit(1);
	}
	continue;
      }else if(*(argv[i]+1)=='D'){
	_nX_=atoi(pc);
	_nY_=_nX_;
	mana = mana | DIMENSION_SET;
	if(_nXmax_<_nX_){
	  fprintf(stderr,"*** Sorry, image dimension cannot be lagaer than %d.\n",_nXmax_);
	  exit(1);
	}else if(_nX_==0){
	  fprintf(stderr,"*** Error: Dimension must not be 0.\n");
	  exit(1);
	}
	continue;
      }else if(*(argv[i]+1)=='N'){
	x=atof(pc);
	Gods_will = Gods_will * pow(10.0,x);
	mana = mana | ENLIGHTENED;
	continue;
      }else if(*(argv[i]+1)=='F'){
	idf=-1;
	n=0;
	while(filterName[n]!=0){
	  if(!strcmp(pc,filterName[n])){
	    idf=n;
	    break;
	  }else{
	    n++;
	  }
	}
	switch(idf){
	case IDF_SHEPP_LOGAN:
	case IDF_RAM_LAK:
	case IDF_RAM_LAK_HANNING:
	case IDF_RAM_LAK_HAMMING:
        case IDF_SHEPP_LOGAN_HANNING:
	case IDF_SHEPP_LOGAN_HAMMING:
	case IDF_NONE:
	  mana = mana | FILTER_SET;
	  continue;
	default:
	  if(0<=idf){
	    fprintf(stderr,
		    "*** Sorry, this filter is not implimented yet: '%s'.\n",
		    filterName[idf]);
	  }else{
	    fprintf(stderr,"*** Error: unrecognized filter name: '%s'\n",pc);
	  }
	  exit(1);
	}
	continue;
      }else{
	pc=argv[i]+1;
	for(pc=argv[i]+1;*pc!=0;pc++){
	  if(*pc=='f'){
	    mana = mana | FORCE_OVERWRITE;
	  }else if(*pc=='v'){
	    mana = mana & (~QUIET);
	    mana = mana | VERBOSE;
	  }else if(*pc=='e'){
	    mana = mana & (~TO_STDOUT);
	    mana = mana | TO_SAME_FILE;
	  }else if(*pc=='I'){
	    mana = mana | TAKE_AREA_SUM;
	  }else if(*pc=='h'){
	    mana = mana | HELP;
	  }else if(*pc=='q'){
	    mana = mana & (~VERBOSE);
	    mana = mana | QUIET;
	  }else if(*pc=='n'){
	    mana = mana | ALLOW_NEGATIVE;
	  }else if(*pc=='c'){
	    mana = mana | KEEP_CASE;
	  }else if(*pc=='p'){
	    mana = mana | FORCE_POSITIVE;
	  }else if(*pc=='a'){
	    mana = mana | DISABLE_SCALE_FCTR;
	  }else if(*pc=='u'){
	    mana = mana | NO_DECAY_CORRECTION;
	  }else if(*pc=='k'){
	    mana = mana | DONOT_RENORM;
	  }else{
            fprintf(stderr,"*** Error: unrecognized option '%s'\n",pc);
	    exit(1);
	  }
	}
	continue;
      }
    }else{
      if(!fnin[0]){
        strcpy(fnin,argv[i]);
        continue;
      }
      if(!fnout[0]){
        strcpy(fnout,argv[i]);
        continue;
      }
    }
    fprintf(stderr,"*** Error: unrecognized argument '%s'\n",argv[i]);
    exit(1);
  }
  if(mana & HELP){
    print_info();
    exit(0);
  }
  if(!fnin[0]){
    fprintf(stderr,"*** Error: sinogram file name missing.\n");
    exit(0);
  }
  if(!fnout[0]&&!(mana & TO_SAME_FILE)){
    fprintf(stderr,"*** Error: output file name missing.\n");
    exit(0);
  }
  if(access(fnin,R_OK)!=0){
    fprintf(stderr,"*** Error: cannot access file: '%s'\n",fnin);
    exit(1);
  }
  if(!(mana & TO_SAME_FILE)){
    if(strlen(fnout)==1 && fnout[0]=='-'){
      mana = mana | TO_STDOUT;
    }else{
      if(access(fnout,R_OK|W_OK)==0){
        if(!(mana & FORCE_OVERWRITE)){
          fprintf(stderr,"*** Error: file already exists: '%s'\n(use -f option to overwrite).\n",fnout);
          exit(1);
        }else{
          if(!(mana & QUIET)) fprintf(stderr,
				     "!!! Warnning: file '%s' will be overwritten.\n",
				     fnout);
        }
      }
    }
  }else if(fnout[0]){
    fprintf(stderr,
	    "*** Error: -e option and output file name cannot be given at the same time.\n");
    exit(1);
  }
  if(mana & VERBOSE) fprintf(stderr,"<< Ground control to USS longshot IV, mission code 0x%x, green?\n",mana);
  if(!(mana & QUIET)){
    if(mana & FILTER_SET)    fprintf(stderr,"(++) Filter No. = %d (%s)\n",idf,filterName[idf]);
    if(mana & ZOOM_SET)      fprintf(stderr,"(++) Zoom = %g\n",zoom);
    if(mana & ROT_SET)       fprintf(stderr,"(++) Rot (deg) = %g\n",rot);
    if(mana & DX_SET)        fprintf(stderr,"(++) dX (cm) = %g\n",DX);
    if(mana & DY_SET)        fprintf(stderr,"(++) dY (cm) = %g\n",DY);
    if(mana & IDX_SET)       fprintf(stderr,"(++) dX (pixel) = %g\n",iDX);
    if(mana & IDY_SET)       fprintf(stderr,"(++) dY (pixel) = %g\n",iDY);
    if(mana & GWcm_SET)      fprintf(stderr,"(++) W (cm) = %g\n",gw);
    if(mana & GWpx_SET)      fprintf(stderr,"(++) W (pixel) = %g\n",gw);
    if(mana & DIMENSION_SET) fprintf(stderr,"(++) Dim (pixel) = %dx%d\n",_nX_,_nY_);
    if(mana & ENLIGHTENED)   fprintf(stderr,"(++) Extra normalization factor = %g\n",Gods_will);
  }
  if(mana & TO_STDOUT) strcpy(fnout,"stdout");
  rot=rot/180.0*PI;
  filter_code=idf;
  if((mana & TAKE_AREA_SUM)){i2=4;}else{i2=1;}
  t0=time(NULL);

  /* USS longshot */
  emi = mana;
 emi:
  init_const();
  if(emi & VERBOSE) fprintf(stderr,
			     ">> emi to ground control, MC 0x%x green, commencing countdown...\n",emi);
  // Open input file
  if(!(mana & QUIET)) fprintf(stderr,"Opening input file: '%s'\n",fnin);
  fin=fopen(fnin,"r");
  if(fin==NULL){
    fprintf(stderr,"*** Error: cannnot open file: '%s'\n",fnin);
    exit(1);
  }
  // Get size of input file.
  fseek(fin,0L,SEEK_END);
  fsin=ftell(fin);
  rewind(fin);
  if(fsin==0){
    fprintf(stderr,"*** Error: Input file is empty: '%s'\n",fnin);
    exit(1);
  }
  if(mana & VERBOSE) fprintf(stderr,"Input file: '%s' %d bytes\n",fnin,fsin);
  // Read main header
  if(mana & VERBOSE) fprintf(stderr,"Reading main header...\n");
  read_header(HID_MAIN,fin);
  if(headers[HID_MAIN][54]!=1){
    fprintf(stderr,"*** Error: File_type = %d is not a scan file.\n",headers[HID_MAIN][54]);
    exit(1);
  }
  read_header_params(HID_MAIN);  
  if(mana & VERBOSE) echo_header_params(HID_MAIN);
  total_number_of_images=guess_number_of_images();
  i=guess_filesize(FID_SCAN);
  if(i<fsin){
    emi = emi | TOO_MANY_SINOGRAM;
  }else if(i>fsin){
    emi = emi | TOO_FEW_SINOGRAM;
  }
  if(emi & (TOO_MANY_SINOGRAM|TOO_FEW_SINOGRAM)){
    if(!(mana & QUIET)) fprintf(stderr,"!!! Warning: total number of sinograms is uncertain.\n!!! Inconsistency between real filesize and expected filesize = %d has been detected;\n!!! if there should be %d sinograms as the main header assures.\n",i,total_number_of_images);
    strcpy(str_total,"??");
  }else{
    sprintf(str_total,"%d",total_number_of_images);
  }
  // Guess output file size just in passing.
  fsout=guess_filesize(FID_IMAGE);
  if(mana & VERBOSE) fprintf(stderr,"Output filesize will be %d bytes\n",fsout);
  // Open output file
  if(mana & TO_SAME_FILE){
    p=(long)&original_file_name[0];
    f2cstr(p,strlen(original_file_name));
    strcpy(fnout,original_file_name);
    p=(long)&fnout[0];
    i=at_exp(p,strlen(fnout));
    pc=&fnout[0];
    pc=pc+i;
    memcpy(pc,".IMG",4);
    pc=pc+5;*pc=0;
    if(!(mana & KEEP_CASE)){
      p=(long)&fnout[0];
      lower(p,strlen(fnout));
    }
    if(access(fnout,R_OK|W_OK)==0){
      if(!(mana & FORCE_OVERWRITE)){
        fprintf(stderr,"*** Error: file already exists: '%s'\n(use -f option to overwrite).\n",fnout);
        exit(1);
      }else{
        if(!(mana & QUIET)) fprintf(stderr,
				 "!!! Warnning: file '%s' will be overwritten.\n",
				 fnout);
      }
    }
  }
  if(!(mana & QUIET)) fprintf(stderr,"Opening output file: '%s'\n",fnout);
  if(mana & TO_STDOUT){
    fout=stdout;
  }else{
    fout=fopen(fnout,"w");
    if(fout==NULL){
      fprintf(stderr,"*** Error: cannot open output file: '%s'\n",fnout);
      exit(1);
    }
  }
  // Set main header for image file
  set_img_header(HID_MAIN);
  // Write main header for image file
  write_header(HID_MAIN,fout);
  // Set some params for FBP calculation
  a=(double)_nx_/(double)_nX_/zoom; // Sample distance to pixel width ratio.
  d=PI/(double)_nt_;                // projection direction difference.
  slice_width=0.675;
  recon_scale=zoom;
  x0=(double)(_jX0_)-iDX-1.0; // origin locates at the center of pixel.
  y0=(double)(_jY0_)-iDY-1.0;
  if(mana & TAKE_AREA_SUM){
    dx[0]=hm*a-x0*a-DX/a;dy[0]=hp*a-y0*a-DY/a;
    dx[1]=hp*a-x0*a-DX/a;dy[1]=hp*a-y0*a-DY/a;
    dx[2]=hm*a-x0*a-DX/a;dy[2]=hm*a-y0*a-DY/a;
    dx[3]=hp*a-x0*a-DX/a;dy[3]=hm*a-y0*a-DY/a;
  }else{
    dx[0]=-x0*a-DX/a;dy[0]=-y0*a-DY/a;
  } 
  // Main loop for each sinogram in a file.
  i_bed_position=0;
  i_planes=0;
  i_frames=1;
  i_record=0;
  i_matrix=0;
  forward_block=2;
  backward_block=0;
  number_of_matrix_created=0;
  number_of_matrix_available=0;
  if(emi & VERBOSE){
    fprintf(stderr,">> This is emi to ground control, EC=0x%x, green, go.\n",emi);
  }
    for(n=1;;n++){
    t1=time(NULL);
    // Read matrix record, if necessary
    if(number_of_matrix_created<=i_matrix) i_matrix=0;
    if(i_matrix==0){
      if(last_record) break;
      i_record++;
      if(mana & VERBOSE) fprintf(stderr,"Reading matrix record No. %d...\n",i_record);
      fread(&headers[HID_MATREC][0],_BS_,1,fin);
      read_header_params(HID_MATREC);
      if(mana & VERBOSE) echo_header_params(HID_MATREC);
      // Set matrix record and write it.
      set_img_header(HID_MATREC);
      fwrite(&headers[HID_MATREC][0],_BS_,1,fout);
    }
    i_matrix++;
    // Get FR BD PL from matrix record
    read_someheader_params(i_matrix);
    // Let's begin.
    if(mana & VERBOSE) fprintf(stderr,">>>> ");
    if(!(mana & QUIET)) fprintf(stderr,"Reconstructing image %2d [FR %2d BD %2d PL %2d] of %s... ",
				n,i_frames,i_bed_position,i_planes,str_total);
    if(mana & VERBOSE) fprintf(stderr,"\n");
    // Read sinogram data and its subheader.
    if(scn2mtv(fin,_nx_,_nt_,_BS_*3,_BS_,n)!=0){
      fprintf(stderr,"*** Error: unexpected EOF detected during read of sinogram\n*** No. %2d [FR %2d BD %2d PL %2d] at record %d\n",
				n,i_frames,i_bed_position,i_planes,i_record);
      fclose(fin);
      if(fout!=stdout) fclose(fout);
      exit(1);
    };
    read_header_params(HID_SUB);
    pixel_size=a*sample_distance;
    Delta=sample_distance;
    if(mana & VERBOSE) fprintf(stderr,"Pixel size = %g cm\n",pixel_size);
    if(mana & VERBOSE) echo_header_params(HID_SUB);
    if(frame_duration==0){
      if(!(mana & VERBOSE)) fprintf(stderr,"\n");
      if(!(mana & QUIET)) fprintf(stderr,"!!! Warnning: Frame_duration = 0 is detected.\n!!! Forced to be 1 and decay correction disabled.\n");
      frame_duration=1;
      mana = mana | NO_DECAY_CORRECTION;
    }
    if(scale_factor==0.0){
      if(!(mana & VERBOSE)) fprintf(stderr,"\n");
      if(!(mana & QUIET)) fprintf(stderr,"!!! Warnning: Scale_factor = 0 is detected.\n!!! Forced to be 1.0\n");
      scale_factor=1.0;
    }
    // Calculate decay correction factor.
    if(!(mana & NO_DECAY_CORRECTION)){
      decay_correction_factor=get_decay_correction_factor();
    }else{
      decay_correction_factor=1.0;
    }
    if(mana & VERBOSE) fprintf(stderr,"Decay correction factor = %g\n",decay_correction_factor);
    // A constant factor appliled to pixel value
    if((mana & DISABLE_SCALE_FCTR)) scale_factor=1.0;
    norm=Delta
        /(2.0*(double)_nt_)
        *(double)decay_correction_factor
        /((double)(frame_duration)*1.0e-3)
        /slice_width
        *pow(pixel_size,2.0)
        *Gods_will
        *scale_factor;
    if(!(mana & DONOT_RENORM)) norm=norm*pow(zoom,2.0)*pow((double)_nX_,2.0);
    if(mana & TAKE_AREA_SUM) norm=norm/4.0;
    // Gaussian smooting.
    if(mana & (GWpx_SET|GWcm_SET)){
      g_width=gw;
      if(mana & GWpx_SET) g_width=g_width*Delta;
      if(mana & VERBOSE) fprintf(stderr,"Gaussian smoothing width = %g cm\n",g_width);
      smtyi();
    }
    // Set Fourier transform of convolution matix.
    if(first) setfx(idf);
    first=0;
    setfq();
    // Clear image matrix
    clelmd();
    // Loop for 4-point integraion, if necessary.
    n_in=0;
    n_out=0;
    for(i=0;i<i2;i++){
      for(jY=0;jY<_nY_;jY++){
	ky=a*(double)jY+dy[i];
	for(jX=0;jX<_nX_;jX++){
	  kx=a*(double)jX+dx[i];
	  s=0.0;
          // summation with respect to projection direction
          for(it=0;it<_nt_;it++){ 
	    x=(double)it*d;
	    r=kx*cos(x+rot)+ky*sin(x+rot);
            // Decomposition of projected coordinate into what-would-u-call-it-um.
	    m0=floor(r);
	    f0=r-(double)m0;
	    mm=m0+_ix0_*_BFS_;
	    n_in++;
	    if(mm<2) continue;
	    if(_nx_*_BFS_-1<=mm) continue;
	    // 5-point Lagrange interpolation of what-did-u-call-it.
	    f0sqr=pow(f0,2.0);
	    s=s+(f0sqr-1.0)*f0*(f0-2.0)/24.0*fq[it][mm-2]
	       -(f0sqr-4.0)*f0*(f0-1.0)/6.0*fq[it][mm-1]
               +(f0sqr-1.0)*(f0sqr-4.0)/4.0*fq[it][mm]
	       -(f0sqr-4.0)*f0*(f0+1.0)/6.0*fq[it][mm+1]
	       +(f0sqr-1.0)*f0*(f0+2.0)/24.0*fq[it][mm+2];
	    n_out++;
	  }
          lmd_j[jY][jX]=lmd_j[jY][jX]+s*norm;
        }
      }
    }
    if(n_in!=n_out){
      if(!(mana & VERBOSE)) fprintf(stderr,"\n");
      if(!(mana & QUIET)) fprintf(stderr,"!!! Warnning: Out-of-index exception has been detected %d times out of %d.\n",n_in-n_out,n_in);
    }
    // float to integer conversion of pixel value.
    quant_scale=f2ilmd(&min,&max,&image_min,&image_max,(mana & FORCE_POSITIVE));
    if(mana & VERBOSE) fprintf(stderr,
			     "Quantification: %g [%g,%g] -> [%d,%d]\n",
			     quant_scale,min,max,image_min,image_max);
    if(!(mana & ALLOW_NEGATIVE)) image_min = 0;
    // get elapsed time in reconstruction.
    recon_duration=(int)(time(NULL)-t1);
    if(!(mana & QUIET)){
      fprintf(stderr,"elapsed: ");
      if(0<recon_duration){
	  fprintf(stderr,"%d",recon_duration);
      }else{
	fprintf(stderr,"1>");
      }
      fprintf(stderr," sec.\n");
    }
    // write image and subheader
    set_img_header(HID_SUB);
    write_data_seq(fout,n);
  }
  fclose(fin);
  if(!(mana & TO_STDOUT)){
    if(!(mana & QUIET)){
      fseek(fout,0L,SEEK_END);
      fsout=ftell(fout);
      fprintf(stderr,"Output file: '%s', %d bytes written.\n",fnout,fsout);
    }
  }
  fclose(fout); 
  if(mana & VERBOSE) fprintf(stderr,"thank you.\n");
  exit(0);
}
